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The present work considers a convection-driven T-shaped standing-wave thermoacoustic system. To gain 
insights on the conversion process of heat to sound and to study the nonlinear coupling between 
unsteady heat release and acoustic disturbances, thermodynamic analysis, numerical and experimental 
investigations are conducted. Three parameters are examined: (1) the inlet flow velocity, (2) heater tem¬ 
perature and (3) heat source location. Their effects on triggering limit cycle oscillations are first investi¬ 
gated in 2D numerical model. As each of the parameters is varied, the head-driven acoustic signature is 
found to change. The main nonlinearity is identified in the heat fluxes. To characterize the transient 
(growing) behavior of the pressure fluctuation, the thermoacoustic mode growth rate is defined and cal¬ 
culated. It is found that the growth rate decreases first and then ‘saturates’. Similar behavior is observed 
by examining the slope of Rayleigh index. Furthermore, the overall efficiency of converting the input 
thermal energy into acoustical energy is defined and calculated. It is found that the energy conversion 
efficiency can be increased by increasing the inlet flow velocity. To validate our numerical findings, a 
cylindrical T-shaped duct made of quartz-glass with a metal gauze attaching on top of a Bunsen burner 
is designed and tested. Supercritical bifurcation is observed. And the experimental measurements show a 
good agreement with the numerical results in terms of mode frequency, mode shape, sound pressure 
level and Hopf bifurcation behavior. 

© 2013 Elsevier Ltd. All rights reserved. 


1. Introduction 

The heat-to-sound or visa versa conversion mechanism (also 
known as thermoacoustics) has fascinated scientists for many 
years now, due to its practical energy applications [1], Standing- 
wave engines is one of the simplest thermoacoustic system, which 
operate with a standing acoustic wave in a resonance tube [2,3], 
The acoustic waves causes pressure oscillations as well as gas 
displacement [4], One of the typical types of resonance tubes is 
Rijke-type combustor [5,6], It is a straight vertical tube with a heat 
source enclosed in its lower half. It provides an elementary 
example of convection-driven thermoacoustic oscillations [2], 
which result from a feedback loop between the heat transferred 
to the fluid and the acoustics propagating in the tube. If the heat 
input is at times of high pressure, a self-amplification of acoustic 
fluctuations may lead to a limit cycle [7,8], The limit cycle is mainly 
due to the nonlinear effects, which cannot be predicted by linear 
theory [5], 

To elucidate or optimize the conversion process in Rijke-type 
thermoacoustic systems, extensive investigations [8-13] have 
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been conducted. The pressure power spectrum of a Rijke combus¬ 
tor [12] is experimentally measured to study the nonlinear 
coupling between heat release rate and the heat-driven sound. 
The coupling process is shown to play an important role in trigger¬ 
ing thermoacoustic oscillations. The transition to Rijke-type ther¬ 
moacoustic oscillations is experimentally investigated [8], To 
predict the stability behaviors of a Rijke tube, theoretical modal 
stability analysis [9] is conducted using Green’s function approach. 
It shows that the thermoacoustic frequency is oscillating with 
increased ‘heater’ location, which is experimentally validated 
[13]. In addition, numerical simulations of Rijke-type thermoacou¬ 
stic oscillations are performed by using FLOW3D [11] or FLUENT 
[10], The oscillations result from the proper ‘phasing’ between 
the heat transferred to the surrounding fluid and the acoustics of 
the resonator. 

Different designs and various forms of input energy in thermoa¬ 
coustic systems have been tested[14-22] to elucidate the heat-dri¬ 
ven large-amplitude or high-frequency sound waves [14,15,23,24]. 
By using an adjustable-power electrical heater, a miniature 
thermoacoustic Stirling engine is designed and tested [16]. 
Thermoacoustic converters powered by engine waste heat [17] or 
solar power [18] have also been developed. Recent developments 
and application of thermoacoustic system for energy harvesting 
are experimentally investigated[3,6]. It is shown that the 
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Nomenclature 


A bottom branch cross-sectional area, m 2 

A t , A r cross-sectional areas of left and right branches, m 2 

Af flame surface area, m 2 

As heater surface area, m 2 

A 0 initial amplitude, m/s 

23 mechanical work done by the fluid supporting a normal 

stress,J 

c speed of sound in air, m/s 

c p , c v heat capacity ratio at constant pressure and volume, kj/ 
kg K 

D the diameter of the bottom stem, m 

£ a generated sound power, J 

(£ q acoustical energy, J 

(£ e exchange rate from heat to sound, J 

Q exponential complex function 

h convection coefficient, W/m 2 I< 

Li, L r left and right branch length, m 

Slf the thermal length scale in the fluid, m 

the thermal length scale in the heating band, m 
M w molecular weight, kg/mol 

n the unit vector of the direction normal to the enclosed 

surface 

p instantaneous pressure, Pa 

Pi, p 0 inlet and outlet pressure, Pa 

q heat flux, W/m 2 

Q heat production per unit volume J/m 3 

Q unsteady heat release rate, kj/s 

Q t input heat power, Watt 

R gas constant, J/mol K 

TZ radius of the tube, m 

9?u universal gas constant, kj/kmol I< 

\>! a normalized Rayleigh index 

S entropy per unit mole, J/mol K 

s specific entropy, J/kg K 

S u flame speed, m/s 

S viscous stress, Pa 

T Temperature, K 


T oscillation period, s 

u(x, t) velocity at a given point x, m/s 
Vi inlet velocity in i direction, m/s 

2B thermal and viscous damping, J 

x coordinates of a given point 

x axial position, m 

Xf heat source location along x 2 , m 

y horizontal axial position, m 

Greek symbols 

a pressure fluctuation amplitude slope, Pa/s 

(5 Kronecker delta function 

Aq r the heat of reaction per unit mass of the mixture, J/kg 

AT S the surface temperature difference, K 

A Tf the temperature difference between the fluid and the 
heater surface, I< 

y the ratio of specific heat 

X conduction coefficient, W/m K 

/( viscosity, kg/m s 

co oscillation frequency, rad/s 

p air density, kg/m 3 

pm the fuel-air mixture density, kg/m 3 

t the viscous stress tensor, N/m 

*P energy conversion efficiency 

a thermoacoustic oscillation growth rate, s~ ] 

Subscript 

s heating band surface 

a, acoustic 

o branch outlet 

Superscript 

mean value 
time derivative 
' fluctuation part 


convection-driven Rijke-Zhao system can harvest 60% more power 
than the conventional conduction-driven standing-wave system 
[3], The standing-wave engine energy conversion efficiency could 
be increased dramatically by using high-pressure gas, such as 
13.8 bar helium [25], However, such engines work less efficient 
than traveling-wave (Stirling) engines. It is experimentally demon¬ 
strated [26] that a traveling-wave engine can achieve 49% of the 
Carnot efficiency. It is also shown [27] that as high as 18% ther- 
mal-to-electric conversion efficiency is achievable by replacing 
piezoelectric generator with electrodynamic alternator. 

There are current interests in more complex geometric configu¬ 
ration of convection-driven thermoacoustic systems. However, few 
designs and investigations are reported. In addition, to optimize 
the conversion from heat to sound, the dynamic interaction be¬ 
tween unsteady heat release [1,2] and acoustic waves needs to 
be better understood. Thus there is a need to identify the roles of 
the heat source and flow conditions, and to investigate the effect 
of its position and the heat fluxes in triggering thermoacoustic 
oscillations, especially in a complex geometric system. Lack of 
these investigations in the literature partially motivated the pres¬ 
ent work. 

In this work, a convection-driven T-shaped thermaocoustic sys¬ 
tem is considered. Both experimental and numerical simulations 
are conducted. To gain insight on the heat-to-sound conversion 


mechanism and the nonlinearity, a general heat-driven acoustic 
wave equation is derived by performing thermodynamic analysis. 
This is described in Section 2. The unsteady Navier-Stokes equa¬ 
tions are iteratively solved by applying a finite volume method 
on spatially discretized mesh of the T-shaped system, as described 
in Section 3. The effects of the inlet flow velocity, heat source 
location and heater temperature on triggering thermoacoustic 
oscillations are then examined. This is described in Section 4. The 
thermoacoustic oscillation growth rate, Rayleigh index, and the 
heat-to-sound energy conversion efficiency are defined and esti¬ 
mated. Finally, experiments are conducted to validate our numer¬ 
ical findings. This is described in Section 5. Comparison is then 
made between the numerical and experimental results. 


2. Thermodynamic analysis 

2.1. Heat-driven acoustic model 

To shed light on thermoacoustic energy conversion, a general 
heat-driven acoustic wave equation is derived by conducting a 
generalized thermodynamic analysis [28]. A T-shaped thermoa¬ 
coustic system is considered, as shown in Fig. 1(a). It consists of 
using multiple concentric constant-temperature rings as a heat 
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(a) 



(b) 

Fig. 1. Schematic of a convection-driven T-shaped standing-wave thermoacoustic system. 


source placed in the bottom stem of the tube. As the air is heated 
up, a buoyancy-driven convection flow with a low Mach number 
M 2 <c 1 is generated. Thus the fluid considered is assumed to sat¬ 
isfy the perfect gas laws. The state equation of perfect gas [29] is 
given as 


R = rt = 3lt = - 

P Mw y 


( 1 ) 


dp _ 1 dp p(y - 1) ds _ 1 dp (y-1) 

dt c 2 dt c 2 dt c 2 dt c 2 ^ [ ’ 

where Q= pTds/dt is the heat addition per unit volume. And it con¬ 
sist of a mean and fluctuating part as Q = Q + Q'. 

The mass and momentum equations over the control volume, 
when linearized with respect to the perturbation quantities, can 
be written as 


where p is the pressure, p density, R is gas constant and T is the 
temperature, 5R„ = 8.3143 kj/kmol I< is the universal constant, 
y = 1.4 is the ratio of specific heat, c is the sound speed, and M w is 
the molecular weight. The pressure, velocity, density and tempera¬ 
ture of the flow at a given location x = (xi, x 2 , x 3 ) = (x, y, z) and time 
consist of a mean part denoted by an overbar, and a fluctuating one 
denoted by a prime as given as 

p(x,t)=p(x) + p'(x,t) u(x, t) = u(x) + u'(x, t) (2) 

p(x,t) = p(x) + p'(x,t) T(x, t) = T(x) + T'(x, t) (3) 

where velocity u(x, t) = (iq, u 2 , u 3 ). 

If the internal energy of the fluid mixture by e, then the first law 
of thermodynamics considering energy conservation can be de¬ 
scribed by 

de pd(p-!) _ ds 

dt + dt dt w 


Now consider a control volume enclosing the heating rings, the den¬ 
sity p varies through changes in both the pressure p and specific en¬ 
tropy s. Thus 


dp 

dt 


dp dp 
dp s dt + 


dp ds 
ds p de 


1 dp dp ds 
c 2 dt + ds p dt 


For a perfect gas p = plRT= yp/c 2 , dpjds j p 
and so Eq. (5) can be rewritten as 


(5) 

-pT(y-l)/c 2 = -p/c p 


dp' _ dll' 
dt 1 dX[ 


(7) 

( 8 ) 


where S denotes the ‘friction'/damping (e.g. viscous stress) effect. 
By substituting linearized Eq. (6) into Eq. (7), multiplying Eq. (7) 
with c 2 p', Eq. (8) with u' and addition, an energy conservation equa¬ 
tion is obtained namely, 


d_ 

dt 


1 c 2 p' : 


2^ +2 p j 


d , , s ,d<5 
+ -(p«5, J U i )-U i - 


‘dx < 


(7-1) 

yp 


p'Q! 


(9) 


Integrating Eq. (9) over the entire system leads to 


d_ 

dt 


1 c 2 p' 2 
2~p~ 


+ 2 PU ' 


dx, = 


where 


( '’ _ p'Q'dXj - S - an (10) 

yp 



(ii) 


It can be seen from Eq. (10) that the term on the left hand side 
describes the classical acoustic disturbance energy (£ a as defined 
as 


(£a = 


T c 2 p' 2 

2 ~W 


2 pU ‘ 


dXi 


( 12 ) 
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The first term on right hand side of Eq. (12) describe the potential 
energy Gj. The other term describes the kinetic energy i.e. 
®a = Puf/2. 

The first term on the right-hand-side of Eq. (10) as defined as (£ e 




( 7 - 1 ) 

yp 


p'Q'dx, 


(13) 


describes the energy exchange between sound and head input. The 
second term S describes the mechanical work done by the fluid 
supporting an externally induced normal stress p' and moving with 
velocity u'. The last term 2B is associated with thermal and viscous 
dissipation.Substituting Eqs. (12) and (13) into Eq. (10) and inte¬ 
grating it over one period leads to 

/ t+T nt+T 

& e dt = & a \ c t +T + J [S + 2B]dt (14) 


Eq. (14) reveals that the input thermal energy is partially converted 
into acoustical energy, partially dissipated by thermal and viscous 
effect, and partially used to produce mechanical work. When there 
is no heat addition, Q = 0 (isentropic), and the working fluid is qui¬ 
escent (stagnant and uniform) and frictionless, Eq. (10) can be sim¬ 
plified into the classical acoustic energy balance equation: 

^ = -» = -||ncW (15) 


where I = p'u' is the acoustic intensity. A is the surface enclosing the 
control volume with outer normal n. Here the Gauss theorem is 
used to transform f v (V ■ I)dV into the surface integral J A (I ■ n)dA 
And it can be used to calculate the acoustic power generated by 
heat addition. 

To produce thermoacoustic oscillations (heat-driven sound), 
one of the necessary condition as suggested by Eq. (14) is 

(£ e > 0, equivalently (16) 


will be produced, even with large amount of heat being added. This 
is known as ‘Rayleigh criterion’. The other critical condition is that 



(£,,dt > 



[SB + SB]dt 


(17) 


That is to say the heat exchange over one period must be greater 
than the sum of that dissipated due to thermal and viscous effects 
and that used to generate mechanical work. To optimize heat- 
to-sound conversion, it is necessary to maximize the difference be¬ 
tween the terms on the left- and right-hand-side in Eq. (17). This 
can be achieved by maximizing (£ e and minimizing Sffi. 

By neglecting the damping term S, substituting Eq. (6) into Eq. 
(7) differentiating Eq. (7) with respect to time t, taking divergence 
of Eq. (8), an inhomogeneous heat-driven wave equation is ob¬ 
tained namely, 


1 d 2 p' 
c 2 dt 2 




(18) 


where overdot denotes the time derivative. Eq. (18) describes a gen¬ 
eral heat-driven acoustic wave equation. Q'(t) describes the time 
rate of heat addition from a heat source. Typically, a flame or an 
electrical heater is used. The instantaneous heat release rate Q from 
such heat sources [5,29,30] can be estimated as 


f p m S u A f Aq r flame 
\ AdT/dxA s heater 


(19) 


where p m is the fuel-air mixture density, S u is the flame speed, Aj is 
the flame surface area, and Aq r is the heat of reaction per unit mass 
of the mixture, 2 is the thermal conductivity, dT/dx is the spatial 
temperature gradient between the heater and surrounding air, A s 
is the surface area of the heater. 


3. T-shaped thermoacoustic model 


that is to say p' and Q are in phase. If |Zp'Q'| > 90°, i.e. p' and Q' are The T-shaped numerical configuration is generated with struc- 

out of phase, then the energy exchange rate (£ e < 0. And no sound tural mesh in Pointwise [31], as shown in Fig. 1(b). Multiple con- 
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Fig. 2. Temperature and velocity contours of the T-shaped thermoacoustic system. 
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Fig. 3. Variation of pressure and heat flux fluctuations with time, as Xf/L = 0.5, 
T s = 1100 K and v, is set to 5 different values. 


branches open ends are pressure outlets p' 0 = 0 Pa. The heating 
band surface temperature is set to 1100 K, which is closed to the 
heat source temperature used in our experiment described later. 
The heat addition causes a convection flow moving towards the 
outlets. 

The steady conservation equations of mass, momentum and en¬ 
ergy with equation of state for ideal gas are solved iteratively for a 
laminar 2D flow as given as 


dp 9(pUj) _ 
dt dXi 
d Uj dUi 0 

p Dt +pui dx i =d X -y pd * +xij) 


fdT 

dT\ 

dUj 

d 

(. dT\ 

V dt + Ui 

dXj) 

]+P dXj- 

= dXj 

{ dxJ 


+ Tij 


dUj 

dXi 


(20a) 

(20b) 

(20c) 


where u, the velocity component of u(x, t) in ith direction, and t,j is 
the viscous stress tensor as given as 


dUi dUj 


2 p du k 

3 dx k ji 


( 21 ) 


Sjj is the Kronecker delta function. When i=j, <5,j = l. Otherwise 
5 t j = 0. p is the dynamic viscosity, and a is the coefficient of heat 
addition to the fluid particle. Q is the instantaneous heat production 
per unit volume from the constant-temperature heating bands. The 
heat transfer rate via conduction over the heating band can be 
shown as 


centric rings with constant surface temperature are equally spaced 
along y-axis. They have a length of L h = 5 cm and placed at Xf. Near 
the heat sources, finer grid (Hexahedra) is used. The geometric 
dimensions are L t = L r = 500 mm, L = 300 mm, A t = A, = TtW, 2 = 1.6 
xl0~ 3 m 2 , A = 2.0 x 10~ 3 m 2 . Before starting the unsteady simu¬ 
lations to capture the thermoacoustic oscillations generated in 
the T-shaped system, a steady solution is first obtained in terms 
of the temperature and velocity contours, as shown in Fig. 2. The 
boundary conditions are set according to our experimental setup 
as shown later as: the bottom stem (mother tube) open end is pres¬ 
sure inlet with a gauge pressure 0.01 Pa, while the bifurcating 


i=l 


1 dxt 


( 22 ) 


where A s is the total surface area of the heating bands. A S i denotes 
the surface area of ith heating cell. As the surrounding fluid flows, 
convection heat transfer is expected. And it is balanced by the con¬ 
duction heat transfer as given as 

Q = hAs(T f - T s ) = hAsAT f = qA s (23) 

where q is the heat flux. It can be shown that the convection heat 
transfer coefficient h can be determined by using 





Time (s) 


Fig. 4. Variation of acoustic power (a), heat input (b), instantaneous velocity and pressure fluctuation (c) and the overall energy efficiency (d). 
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h = X?Ij^ = r( Pr.Re) (24) 

It is known that the convection coefficient is also a function of Pra- 
ndtl and Reynolds number, as shown in Eq. (24). For given charac¬ 
teristic thermal length scales 2/ in the fluid and 2 S in the heating 
band [11], nondimensional Prandtl Pr, Reynolds number Re, and 
Nusssel Nu can be shown as 


Pr = C -f, Re = f,W . 

k /.( k A Tf 


(25) 


where X is the thermal conductivity. In our work, Pr is approxi¬ 
mately 0.86. Rea; 1800 <2300 at the inlet of the bottom branch. 
Thus the inlet flow is a laminar one. 

To solve these conservation equations, 2nd-order schemes are 
chosen for pressure interpolation and spatial discretization, while 
the SIMPLE algorithm is used for velocity-pressure coupling. To 
take the influence of temperature on the physical properties of 
air into account, the specific heat c p (T) = ( 8U/dT ) p , the viscosity 

fi(T) = /.i Kf (r/T re f) 3/2 (here C is Sutherland’s constant) and 
the thermal conductivity x(T) = c p (T)p(T)/Pr are modeled as 
temperature-dependent[32[. The following values at 300 I< are ap¬ 
plied: c p = 1.006 kj/kg K, n = 1.724 x 1(T 5 kg/m s and X = 0.0242 W/ 
m K. 


4. Numerical results 

To gain insights on the interaction between unsteady heat re¬ 
lease and acoustics, parametric studies are conducted. This include 
examining the effects of the inlet flow velocity, heat source loca¬ 
tion, and heater surface temperature on triggering thermoacoustic 
oscillations. These parameters are chosen, based on the facts: (1) 
varying the inlet velocity results in the changes of the convection 
coefficient h = F(Pr, Re) and so the heat flux q, (2) the heat source 
location play an important role in determining the coupling 
between p' and q' (i.e. the phase Zp'q') and so (E e , (3) the heater 
surface temperature is critical in determining the heat flux 
q oc AT = T s — Tf. Each of the parameters is varied, while all the 
other parameters remain constant. 

4.3. Inlet flow velocity effect 

The effect of the inlet Flow Velocity at the bottom open end on 
the resulting thermoacoustic oscillations is first investigated. The 
inlet flow velocity can be varied by changing the bottom open 
end inlet pressure. Fig. 3 shows the variation of the pressure oscil¬ 
lation and heat fluxes, as the inlet velocity is set to 5 different 
values. It can be seen that when the inlet velocity is Vi ^ 0.04 m/ 
s, the initial disturbances decay. However, as the inlet velocity is 
increased, limit cycles occur. And the larger inlet velocity is, the 
larger the resulting pressure and heat flux oscillation amplitude is. 

In order to gain insight on the heat-to-sound energy conversion, 
the total acoustical energy can be estimated by using Eq. (12). 
Fig. 4(a) shows the variation of the acoustic energy with time, as 
the inlet velocity is set to 0.81 m/s. It can be seen that the acousti¬ 
cal power is increased with time by extracting more thermal 
energy. When a limit cycle occurs at f ^ 1.0 s, the acoustic power 
keeps almost constant. This is due to the balanced periodic energy 
exchange between the heat and sound. 

If we define the energy conversion efficiency X F of converting 
input heat power Q t into output acoustical power £ a as 



Fig. 5. Variation of Nusselt number with normalized time t/T 0 (T 0 = 2nla> 0 ). 


then the fraction of heat input being converted into acoustical en¬ 
ergy is shown in Fig. 4(d), as denoted by the solid line. Here p' de¬ 
notes the pressure fluctuation, A is the cross-sectional area of the 
tube, H s is the heating band height, v is the velocity fluctuation, 
Hi denotes the radius of the ith heating ring, and q' is the heat flux 
fluctuation. It can be seen that about 0.79% of total heat addition is 
converted into acoustical energy, which is much efficient than 
0.28% when the inlet velocity is set to 0.21 m/s, as denoted by the 
dot line. Thus the increased pressure fluctuation amplitude with 
the inlet velocity, as shown in Fig. 3(a) indicates the intensified en¬ 
ergy conversion between the heat and sound. 

Flow reversal [5] would occur when the velocity perturbations 
exceed the mean. This is accompanied by large fluctuations in heat 
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Fig. 6. Variation of pressure and heat flux fluctuations, as T s = 1100K, p, = 0.3 Pa, 
p 0 = 0 Pa, the heater is placed at 4 different locations and the corresponding mean 
flow velocities at the bottom branch inlet are 0.30 m/s, 0.56 m/s, 0.84 m/s and 
0.96 m/s: (a) pressure fluctuations (b) heat flux fluctuations. 
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Fig. 7. (a) Variation of the mode growth rate a, as the heater is placed at 4 locations, 
(b) Rayleigh index between q(t)/q and p(X/, t)/p, as xrfL = 3/6. 

release rate. When the reverse flow occurs, the working fluid is 
heated part of the way to the heater temperature on its first pass 
through the heater, and receives more energy on its second pass. 
The variation of the instantaneous velocity between two neighbor¬ 
ing heating bands is shown in Fig. 4(c). The dash line denotes the 
mean flow speed. It can be seen that the fluctuating part is larger 
than the mean part. Furthermore, the presence of the spurs in 
the velocity indicates the flow reversal and nonlinearity of the 
thermoacoustic system. Similar nonlinearity behavior is observed 
in the thermal power added to the system, as shown in Fig. 4(b). 
However, this is different from the pressure fluctuations, as de¬ 
noted by the dash line. 

l°SiolQ'l 

(W/m 2 ) 

I 

^ 0 

-2 

w 

™ -6 

12 3 4 

togiolP'l 
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Fig. 8. Spectra of the heat flux (a) and the pressure fluctuation (b), as x f IL = 3/6, 
T s = 1100 K, Pi = 0.2 Pa. 




As the inlet flow velocity is changed, Nusselt number is found to 
vary as shown in Fig. 5. It can be seen that with inlet flow velocity 
increased, the mean Nu number is increased. This reveals the 
dominant role of convection over the conduction. Furthermore, 
compared with the steady case (corresponding to no thermoacou¬ 
stic oscillations present) as denoted by the dash line, the presence 
of thermoacoustic oscillations results in significant increase of heat 
transfer, which has the potential for energy efficient drying. 


4.2. Heater location effect 


The effect of the heater location on the resulting thermoacoustic 
oscillations is then investigated. Fig. 6 shows the variation of the 
pressure and heat flux fluctuations, as the heater is placed at 4 
different locations. It can be seen that when X//L = 1/6, initial flow 
disturbances decay. However, as Xf/L is increased, the amplitudes 
of the pressure oscillation and heat flux are increased, and limit cy¬ 
cles eventually occur. This is due to the intensified energy conver¬ 
sion from heat to sound. In addition, the growing of the pressure 
oscillations at 0 ^ t < 1.0 as shown in Fig. 6(a) can be approximated 
by using the real part of a complex exponential function 

G{co, a) = «{Ae iMt+,rt } (27) 


where Ao denotes the initial amplitude at t = 0, m denotes the oscil¬ 
lation mode frequency and a describes the mode growth rate. The 
thick dash line in Fig. 6(a) describes s Ji{26.67e i2 ' r216t+4287t }. Compar¬ 
ing it with the amplification of the pressure oscillations as Xf/L = 4/6 
reveals that the exponential function can approximate the growing 
process very well when t< 0.55 s. However, if the growth rate a 
keeps constant, the approximation at t > 0.55 s does not agree with 
the pressure fluctuation growing. This indicates that the growth 
rate a is varied during the amplification process. It is worth noting 
that the heat flux amplification cannot be simply approximated by 
using an exponential function due to its inherent nonlinearity as 
discuss later. 

The time-varying growth rate of the dominant mode of the 
pressure oscillations is shown in Fig. 7. It can be seen that the 





q(t) (W/m 2 ) P'(t) (Pa) 


Fig. 9. Phase-space plots for the pressure oscillations and heat fluxes (a) phase 
diagram between q(t ) and O denotes T s = 700 K, and o denotes T s = 600 K, (b) 
phase diagram between q(t) and p'(t ) T s = 800 K (c) phase diagram between q(t ) and 
q(t) and (d) phase diagram between p'{t) and p'(t): — o - denotes T s = 900 K, and - 
x - denotes T s = 1000 K. 
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Table 1 

Measured growth rate and frequency of the dominant thermoacoustic mode, as Xj/L = 3/6 and Pi = 0.1 Pa. 
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Heater temperature (K) 

600 

700 

800 

900 

1000 

Mode growth rate 

-17.67-> 0 

—9.00-> 0 

0.73-> 0 

1.33-f 0 

1.86-f 

Mode Frequency (Hz) 

NA 

NA 

214.4 

216.38 

218.82 


growth rate is a = tan a < 0 and the pressure fluctuation decays at 
Xf/L = 1/6, since a. is negative (anti-clockwise). However, when the 
heater is moved to Xf/L > 2/6, the mode growth rate <r = tan a > 0 
and the pressure oscillations are amplified. Furthermore, c(x// 
L = 2/6) ^ cr(xp/L = 3/6) ^ ff(xp/L = 4/6). The larger a mean that 
initial triggering disturbances will grow into a limit cycle more 
quickly. Furthermore, a -> 0 with f -> oo, and the mode growth rate 
is saturated. 

Rayleigh index 9}* is generally used as a critical indicator of the 
coupling between heat and sound [5,7,30]. It describes the energy 
exchange intensity between unsteady heat release and the pres¬ 
sure oscillations. If 9b, > 0, the acoustical energy in the present sys¬ 
tem is increased and tends to drive the unstable process into a 
limit cycle. However, a negative 9U indicates a ‘destruction’ inter¬ 
action between unsteady heat release and pressure oscillations and 
so the thermoacoustic oscillations decay with time. Thus it would 
be interesting to check the Rayleigh index 9i„. Following the work 
[24J, we defined a normalized Rayleigh index 9hi as given as 


« A (t) = JJ <5(x - Xp)dx 


_ p'(m, i k)q'(co, I/O# 

sjf 0 n/to p*(a), W /T 0 27C/[u *)d* 


dt j/ 


(28) 


Fig. 7(b) shows the variation of the Rayleigh index 9i^(t) with time, 
as Xf/L = 3/6. It can be seen that when initial perturbation grows into 
a limit cycle, Rayleigh index 9?* is positive, keeps increasing and 
finally ‘saturated’ at approximately 0.90. This verifies our numerical 
model and explains why self-sustained thermoacoustic oscillations 
occur in the T-shaped system. 

Fig. 8 shows the time-frequency analysis of the heat flux q'(t) 
and the pressure oscillations p'(t). It can be seen that the spectra 
of the heat flux consists of one dominant frequency at approxi¬ 
mately 220 Hz and two harmonics with large amplitude. However, 
the harmonics are negligible in the pressure spectra. This confirms 
that acoustic waves can be approximated by using linearization 
analysis. However, nonlinearity is mainly arising from the unstea¬ 
dy heat release. 


growth rate means the pressure and heat flux fluctuations are 
decaying. And the system is stable with no thermoacoustic oscilla¬ 
tions generated. However, when the growth rate is positive, the 
fluctuations are amplified and finally grow into a limit cycle. The 
frequency of the dominant mode is found to slightly shifted, as 


Fig. 10. Experimental setup of T-shaped thermoacoustic system with a fine metal 
gauze attached on top of a Bunsen burner. 


(a) 




4.3. Heat flux effect 

As shown in Eq. (13), the heat flux q(t) <x(T s - Tf) plays an 
important role in generating thermoacoustic oscillations. To inves¬ 
tigate its effect, the heater temperature T s is set to 5 different 
values and so the mean heat addition q. Fig. 9 illustrates the vari¬ 
ations of the heat flux q(t) and pressure fluctuations p'(t), as the 
heating bands are set to different temperatures. It can be seen that 
the mean heat flux is increased with increased heater temperature. 
Furthermore, the instantaneous heat flux q(t ) and pressure fluctu¬ 
ations p'(t) decay, as the heater temperature is lower than T s - 
^ 700 K. However, as the temperature is above 800 K, as show in 
Fig. 9(b)-(d), the nonlinear thermoacoustic oscillations as indi¬ 
cated by the circle-shaped diagrams are generated. And the ampli¬ 
tudes of the heat flux and the pressure fluctuations are increased 
with the heater temperature. Moreover, the main nonlinearity is 
identified in the rate of heat release, as revealed by comparing 
the phase diagrams of the heat fluxes (pomelo-like) and pressure 
fluctuation (circle-shaped). 

In addition, the growth rate of the dominant thermoacoustic 
mode is calculated and shown in Table 1. The negative mode 






Frequency (Hz) 


Fig. 11. Measured pressure oscillations, as the fuel flow rate is set to 230 ml/s, and 
the heat source is placed at 4 different locations (a) time-domain measurement, 
(b)-(e) frequency spectra at XjjL = 1 /6, XjjL = 2/6, x f /L = 3/6 and Xj/L = 4/6 
respectively. 
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Fig. 12. Temperature contour measurements by using an infrared camera, as the heat source is placed at 4 different locations. 

the heater temperature is set to different values. This might be due involved with 3 microphones. The geometry, dimensions of the 

to the variation of the sound speed c oc Vf. setup and the flow parameters are shown in Appendix. 


5. Description of experiment 


5.2. Experimental results 


5.1. Experimental setup 


In order to validate our numerical findings, a T-shaped tube 
with a premixed laminar flame confined is designed and tested. 
It is made of quartz glass. And it has a mother tube (bottom stem), 
which splits into two daughter tubes (i.e. horizontal branches) 
with equal lengths, see Fig. 10. T-shaped junctions protruding from 
both bifurcating branches are designed for placing microphones. 
To measure pressure fluctuations along the setup, two arrays of 
B&K microphones (Type 4794) are implemented: each array is 
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Fig. 13. Comparison of the 
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and predicted mode shapes. 


5.2.1. Parametric measurements 

Experimental measurements are conducted in two sets, corre¬ 
sponding to our numerical simulations. One set of measurements 
are performed by placing the flame at 4 different locations, as 



p’(t) p’(t) 


Fig. 14. Measured pressure oscillations, as the fuel flow rate is set to 4 different 
values: (a) 150 cc/min, (b) 200 cc/min, (c) 230 cc/min and (d) 260 cc/min. 
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Fig. 15. Temperature contour measurement by using an infrared thermal imaging camera, as the fuel flow rate is set to 4 different values. 


shown in Fig. 11(a) in time domain. Fig. ll(b)-(e) shows the fre¬ 
quency spectra of measured and predicted pressure fluctuations. 
It can be seen from Fig. 11(a) that thermoacoustic oscillations 
could be generated, depending on the flame location. Correspond¬ 
ingly, the mean flow velocity resulting from natural convection at 
the inlet of the bottom stem is increased from 0.04 m/s, 0.10 m/s, 
0.17 m/s to 0.25 m/s as measured by using a thermal anemometer 
TA5, as the flame location is varied from 1/6 to 4/6. When a limit 
cycle is produced, the maximum pressure oscillation can reach 
about 136 dB, as shown in Fig. 11(b). The dominant mode is at 
around 216 Hz and it has multiple harmonics (indicating nonlin¬ 
earity) as shown in the graphs (c)-(e). Comparison of the experi¬ 
mental measurements and the numerical predictions shows good 
qualitative agreement in terms of the sound pressure level (SPL) 
and the mode frequencies. 

The temperature along the T-shaped system is measured by 
using an infrared camera [4,6]. It captures the quartz tube surface 
temperature spectrum. However, it can represent the ‘hot’ status of 
the sound waves inside each branch. The camera is set to operating 
range of 0-500° Celsius (corresponding to 273-773 K). Fig. 12 
shows the measured temperature contour. It can be seen that there 
is a local maximum temperature at the junction section. And the 
temperature is decreased along the axis of the horizontal branch. 
The measured temperature contour (mushroom-shaped) is similar 
as our numerical prediction as shown in the top graph of Fig. 2. 

To gain insight on the thermoacoustic oscillations, the mode 
shapes along the horizontal branches are measured, as shown in 
Fig. 13(a). It can be seen that pressure nodes are located at each 
ends. And the maximum pressure fluctuation occurs at approxi¬ 
mately 0.4 m away from the open ends. Comparing the measured 
mode shape with the numerical predicted one, good agreement 
is obtained. Fig. 13(b) shows the mode shape along the bottom 
stem and the right branch. The total length is corresponding to half 
of wave length of the thermoacoustic oscillations at approximately 
220 Hz. It is worth noting that the mode shapes are measured by 
using the classical two-microphone technique [33]. The sound 
speed is obtained by using the mean value of the local temperature 
measurements. And the mean flow speed is neglected. 



Fig. 16. Hopf bifurcation diagram of measured and numerical predicted thermoa¬ 
coustic oscillations, as the fuel flow rate is set to 300 cc/min. 


The other set of experiments is to investigate the effect of 
heat flux. This is achieved by varying the fuel flow rate. Fig. 14 
shows the measured pressure oscillations, as the fuel flow rate 
is set to 4 different values. It can be seen that when the fuel flow 
rate is lower than 150 cc/min, no limit cycle is generated. This is 
most likely due to the fact that the thermal energy input is not 
larger than that dissipated due to viscous and thermal effects. 
The temperature contour measurement as shown in Fig. 15 is 
similar to those as shown in Fig. 12. However, the local maxi¬ 
mum temperature at the junction section is decreased from 
711.75 K to 685.55 K, as the fuel flow rate is increased. This is 
most likely due to the intensified heat-to-sound conversion (see 
Fig. 16). 











346 


S. Li, D. Zhao/Energy Conversion and Management 75 (2013) 336-347 


5.2.2. Supercritical Hopf bifurcation 

A bifurcation diagram in terms of the pressure fluctuation 
amplitude can be obtained by varying the heater location gradually 
while all the other parameters remain constant. The nonlinear 
behavior around the Hopf bifurcating point xj is characterized by 
a gradually increase of the amplitude as x f > x c f . ForXf < x^, all per¬ 
turbations imposed on the thermoacoustic system tend to decay to 
zero. 

In general, the measured results are in good agreement with the 
numerical predictions, in terms of the mode shapes, mode frequen¬ 
cies and the supercritical Hopf bifurcation diagram. This reveal that 
the numerical model developed can be used to study the onset of 
thermoacoustic oscillations and facilitate the design of convec¬ 
tion-driven standing-wave heat engine systems, which might be 
utilized for energy harvesting or drying. 

6. Discussion and conclusions 

In summary, a nonlinear T-shaped standing-wave thermoacou¬ 
stic system is numerically and experimentally investigated to gain 
insights on the energy conversion process from heat to sound. The 
physical processes of the generation and propagation of sound 
waves by inputting thermal energy is mathematically illustrated 
by performing a generalized thermodynamic analysis. Three key 
parameters are identified and examined: (1) inlet flow velocity, 
(2) the heater location, and (3) heater surface temperature. Their 
effects on triggering limit cycle oscillations are first investigated 
in 2D numerical model by varying one parameter, while all the 
other parameters remain constant. The main nonlinearity is identi¬ 
fied in the heat fluxes. 

To characterize the transient (growing) behavior of the pressure 
fluctuation, the thermoacoustic mode growth rate is defined and 
calculated. It is found that the growth rate decreases first and then 
‘saturates’. Similar behavior is observed in examining the slope of 
Rayleigh index, which is widely used as an onset/triggering indica¬ 
tor to characterize the phasing between unsteady heat release and 
sound. The transient behavior of heat flux cannot be simply de¬ 
scribed by a growth rate due to its nonlinearity. In addition, the 
overall efficiency of converting the input heat power into acoustic 
power is estimated and compared. It is found that the energy con¬ 
version efficiency can be increased by increasing the inlet flow 
velocity. 

To validate our numerical findings, a cylindrical T-shaped duct 
made of quartz-glass with a metal gauze attaching on top of a Bun¬ 
sen burner is designed and tested. To monitor the temperature and 
acoustic fields in the bifurcating branches, an infrared thermal 
imaging camera and two arrays of microphones are used. As a pre¬ 
mixed laminar flame is placed in the bottom branch of the tube, 
self-sustained thermoacoustic oscillations might be generated, 
depending on the 3 identified critical parameters. It is shown that 
Hopf bifurcation is associated with the thermoacoustic system. 


Table 2 

Geometry and flow parameters of the T-shaped thermoacoustic system. 


Geometry 

Value 

Left branch length, L/ (m) 

0.5 

Right branch length, L r (m) 

0.5 

Bottom branch length, L (m) 

0.3 

Cross-sectional area of left branch, A/ (m 2 ) 

0.0016 

Cross-sectional area of right branch, A r (m 2 ) 

0.0016 

Cross-sectional area of bottom branch, A (m 2 ) 

0.002 

Inlet temperature, T u (K) 

298 

Inlet and outlet mean pressure, p u = p d (Pa) 

1.01325 x 10 5 

Fuel, 

Propane 

Lower Heating value, Q (MJ/kg) 

46.35 


Comparison is then made between the experimental measurement 
and numerical predictions. Good quantitative agreement is ob¬ 
tained in terms of the mode shape, mode frequency, sound pres¬ 
sure level (SPL) and the behavior of supercritical bifurcation. 
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Appendix A 

The geometry, dimensions of the setup and the flow parameters 
are shown in Table 2. 
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